#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdio>
#include <string>
#include <iomanip>

#include "B2_function.h"
#include "shekel_function.h"
#include "griewank_function.h"
#include "rosenbrock_function.h"
#include "zakharov_function.h"

#include "continuous_tabusearch.h"
#include "Time.h"
#include "nr_c/nr.h"
#include "gsl_stuff/gsl_function_setup.h"
//#include "wrapperFunctions.h"

#define PI   3.1415926535897932385

using namespace std;

void run_GSL_CGM(int argc, char *argv[ ], int seed, double ftol, Function *func, int nruns, int numdim, string& directory)
{
  ofstream fpParams("./output/CGM_params.used");
  fpParams << "Seed               " << seed << endl;
  fpParams << "Tolerance          " << ftol << endl;
  fpParams << "N_runs             " << nruns << endl;
  fpParams << "N                  " << numdim << endl;

  fpParams.close();

  initiate_ran1(seed);      // initiates the random number

  Time timer;               // to measure calculation times
  cout << "prepare timer " << timer;

  double *sol;              // the solutions we will follow
  //double *deriv;
  double fTol = ftol;
  int nRuns = nruns;
  int numDim = numdim;


  ////////////////////////////////////////////////////////////////
  // Open files CGM_Solution.out and CGM_initial_values.out,
  // and prepare them for output. CGM_Solution.out holds the
  // best function value after each run both without local
  // minimisation and with local minimisation. CGM_initial_values
  // holds the initial coordinates (debugging).
  ////////////////////////////////////////////////////////////////

  ofstream outEndCGM_Solution( string(directory + "CGM_Solution.out").c_str());

  outEndCGM_Solution << "Dimension: " << numDim << endl
    << "Smallest function value after each run\n\n"
    << "      Func Val Bef     Func Val After      Function Evaluations      Distance       Run Number" << endl
    << endl;

  //ofstream finalDeriv( string(directory + "CGM_EndDerivs.deriv").c_str());
  //ofstream finalCoord( string(directory + "CGM_EndCoord.coord").c_str());

  ofstream outIniVal("output/CGM_initial_values.out");

  Function *function;       // Creates an instance of the function
  function = func;          // Sets the function = to the function passed as an argument

  newA(sol, numDim);
  //newA(deriv, numDim);

  cout << "Time: " << timer;

  /*ftolA[1]=1e-50;
  ftolA[2]=1e-40;
  ftolA[3]=1e-30;
  ftolA[4]=1e-20;
  ftolA[5]=1e-10;
  ftolA[6]=1e-5;
  ftolA[7]=1e-4;
  ftolA[8]=5e-4;
  ftolA[9]=1e-3;
  ftolA[10]=1e-2;
  ftolA[11]=1e-1;*/

  //ftol = ftolA[1];


  //////////////////////// Prepare GSL stuff

  size_t iter = 0;
  int status;
  
  const gsl_multimin_fdfminimizer_type *T;
  gsl_multimin_fdfminimizer *gslMinimizer;

  gsl_multimin_function_fdf gslFunc;
  gslFunc.n = numDim;
  gslFunc.f = &gsl_f;
  gslFunc.df = &gsl_df;
  gslFunc.fdf = &gsl_fdf;
  gslFunc.params = function;

  gsl_vector* gslVector;
  gslVector = gsl_vector_alloc(numDim);
  
  T = gsl_multimin_fdfminimizer_conjugate_fr;
  gslMinimizer = gsl_multimin_fdfminimizer_alloc(T, numDim);

  gsl_multimin_fdfminimizer_set(gslMinimizer, &gslFunc, gslVector, 0.01, 1e-4);

  //////////////////////// end prepare GSL stuff


  for(int i=1; i<=nRuns; i++)  
  {
    cout << "Run No. " << i << endl;

    //ftol = ftolA[i];

    function->nCall = 0;
    // Populate the Input array with a random start position/configuration
    function->createRandomSol(sol);

    //finalDeriv << i << setw(14);

    // Print out the run number and initial
    // starting configuration (coordinates)
    outIniVal << "\nRun No. " << i << endl << endl;
    outIniVal << "Variable | Starting values" << endl;
    for (int k=1; k<=numDim; k++)
    {
      outIniVal << setw(5) << k << setw(18) << sol[k] << endl;
    }


    ///////////////////// Set GSL starting point

    for (int k=1; k<=numDim; k++)
      gsl_vector_set(gslVector, k-1, sol[k]);



    // reset tracker before starting a new CTS run (in run_test.h)
    //ResetTracker(sol, function);            // calculates the function value of sol

    double origValue = function->Value(sol);

    // print out best solution (before local minimisation)
    /*char resultCGM[80];
    char resultCGMpstruc[80];
    
    sprintf_s(resultCGM, 80, string(directory + "result_CGM%d.txt").c_str(), i);
    function->print_sol(sol, origValue, resultCGM, 0, 0);*/

    // print out best solution so far (without local minimisation)
    outEndCGM_Solution << setw(16) << origValue;

    /////////////////// do GSL local minimization

    do 
    {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate(gslMinimizer);

      if (status)
        break;

      status = gsl_multimin_test_gradient (gslMinimizer->gradient, 1e-3);

      if (status == GSL_SUCCESS)
        cout << "\n\nSuccess\n\n";
    }
    while (status == GSL_CONTINUE && iter < 100);



    // after local minimisation lBestFuncValue has been changed
    // this sets lBestFuncValue equal to bestFuncValue which is
    // then outputted
    outEndCGM_Solution << setw(18) << origValue << setw(23) << function->nCall << setw(19)
      << function->distanceFromGlobalMinimum(sol) << setw(14) << i << setw(14) << fTol << endl;
  }

  gsl_multimin_fdfminimizer_free(gslMinimizer);
  gsl_vector_free(gslVector);

  double *bestSol;
  newA(bestSol, numDim);
  char resultBest[80];

  sprintf_s(resultBest, 80, string(directory + "CGM_best_struc.txt").c_str());
  double bestVal = function->BestSoFar(bestSol);
  function->print_sol(bestSol, bestVal, resultBest, 0, 0);

  outEndCGM_Solution << "\nGlobal best value over all runs:\t" << setw(9) << function->BestSoFar() << endl;

  //finalCoord.close();
  outEndCGM_Solution.close();
  outIniVal.close();

  // print final output to screen
  cout << endl << endl << "Output files have been written to " << directory << endl;
  cout << "The total time for the simulation was " << timer << endl << endl;
  cout << "---------------------------" << endl << endl;
}

